Method and apparatus for specifying and visualizing robust tuning of model-based controllers

ABSTRACT

A method includes obtaining information identifying uncertainties associated with multiple parameters of a model for an industrial model-based controller. The method also includes obtaining information identifying multiple tuning parameters for the controller. The method further includes generating a graphical display identifying (i) one or more expected step responses of an industrial process that are based on the tuning parameters of the controller and (ii) an envelope around the one or more expected step responses that is based on the uncertainties associated with the parameters of the model. The parameters could include a process gain, a time constant, and a time delay associated with the model. The uncertainties associated with the parameters of the model could include, for each parameter of the model, an uncertainty expressed in the time domain. The information identifying the tuning parameters could include a settling time and an overshoot associated with the controller.

CROSS-REFERENCE TO RELATED APPLICATION AND PRIORITY CLAIM

This application claims priority under 35 U.S.C. §119(e) to U.S. Provisional Patent Application No. 61/954,912 filed on Mar. 18, 2014. This provisional patent application is hereby incorporated by reference in its entirety.

TECHNICAL FIELD

This disclosure relates generally to industrial process control systems. More specifically, this disclosure relates to a method and apparatus for specifying and visualizing robust tuning of model-based controllers.

BACKGROUND

Model predictive control (MPC) techniques use one or more models to predict the future behavior of an industrial process. Control signals for adjusting the industrial process are then generated based on the predicted behavior. MPC techniques have become widely accepted in various industries, such as the oil and gas, pulp and paper, food processing, and chemical industries.

When tuning an MPC or other model-based process controller for industrial use, it is often necessary or desirable to find tuning parameters that ensure good performance in spite of both (i) process disturbances and (ii) mismatches between a model used by the controller and the actual process. This problem falls into the discipline of “control theory” and the practice known as “robust control.” Standard robust control techniques use a concept known as “unstructured uncertainty,” which generally involves analyzing and specifying performance in the frequency domain. These robust control techniques have been used to successfully tune process controllers in a variety of industries.

SUMMARY

This disclosure provides a method and apparatus for specifying and visualizing robust tuning of model-based controllers.

In a first embodiment, a method includes obtaining information identifying uncertainties associated with multiple parameters of a model for an industrial model-based controller. The method also includes obtaining information identifying multiple tuning parameters for the controller. The method further includes generating a graphical display identifying (i) one or more expected step responses of an industrial process that are based on the tuning parameters of the controller and (ii) an envelope around the one or more expected step responses that is based on the uncertainties associated with the parameters of the model.

In a second embodiment, an apparatus includes at least one memory configured to store (i) information identifying uncertainties associated with multiple parameters of a model for an industrial model-based controller and (ii) information identifying multiple tuning parameters for the controller. The apparatus also includes at least one processing device configured to generate a graphical display that identifies (i) one or more expected step responses of an industrial process that are based on the tuning parameters of the controller and (ii) an envelope around the one or more expected step responses that is based on the uncertainties associated with the parameters of the model.

In a third embodiment, a method includes obtaining information identifying uncertainties associated with multiple parameters of a model for an industrial model-based controller. The method also includes identifying multiple tuning parameters for the controller based on the uncertainties and one or more tuning specifications, where the one or more tuning specifications include constraints on a settling time and an overshoot associated with the controller.

Other technical features may be readily apparent to one skilled in the art from the following figures, descriptions, and claims.

BRIEF DESCRIPTION OF THE DRAWINGS

For a more complete understanding of this disclosure, reference is now made to the following description, taken in conjunction with the accompanying drawings, in which:

FIG. 1 illustrates an example web manufacturing or processing system according to this disclosure;

FIG. 2 illustrates an example internal model control structure employed for machine direction model predictive control (MD-MPC) according to this disclosure; and

FIGS. 3 through 12 illustrate details of an example technique for specifying and visualizing robust tuning of model-based controllers according to this disclosure.

DETAILED DESCRIPTION

FIGS. 1 through 12, discussed below, and the various embodiments used to describe the principles of the present invention in this patent document are by way of illustration only and should not be construed in any way to limit the scope of the invention. Those skilled in the art will understand that the principles of the invention may be implemented in any type of suitably arranged device or system.

As noted above, standard robust control techniques have been used to successfully tune model predictive control (MPC) and other model-based controllers in a variety of industries. However, standard robust control techniques often rely upon highly-trained and highly-knowledgeable personnel. This often makes it more difficult and expensive to tune model-based controllers.

This disclosure provides techniques that support robust tuning of MPC and other model-based controllers. Among other things, these techniques can be used by non-expert users or other users to specify model uncertainty and controller performance (i) in the time domain and (ii) with reference to process step responses. While unstructured uncertainty is not an intuitive concept, non-expert or other users can specify a range of uncertainty for the values of simple model parameters. This allows the users to specify robust control designs using concepts that are intuitive and easy to understand.

Depending on the implementation, the techniques disclosed in this patent document allow for the robust tuning of model-based controllers. Example features of the techniques can include:

(1) the specification of model uncertainty in terms of simple model parameters;

(2) the specification of controller performance using simple time domain concepts (such as settling time and overshoot);

(3) the visualization of robust performance through step and disturbance response plots with envelopes of possible responses given the uncertainty of a model;

(4) an algorithm (such as a MATLAB algorithm) that takes user-friendly specifications and returns appropriate tuning parameters; and

-   -   (5) an algorithm (such as a MATLAB algorithm) that provides         robust performance envelope visualization.

In some embodiments, a user interface can be provided that allows users to enter model uncertainty specifications and performance specifications. The user interface can also allow users to view resulting tuning parameters and visualize the resulting controller performance.

Note that in the following description, an example of this functionality is given with respect to use with a controller in a paper manufacturing system. However, this disclosure is not limited to use with controllers in paper manufacturing systems. The techniques disclosed in this patent document can be used with any suitable model-based controller that is used to control any aspect(s) of a process.

FIG. 1 illustrates an example web manufacturing or processing system 100 according to this disclosure. As shown in FIG. 1, the system 100 includes a paper machine 102, a controller 104, and a network 106. The paper machine 102 includes various components used to produce a paper product, namely a paper web 108 that is collected at a reel 110. The controller 104 monitors and controls the operation of the paper machine 102, which may help to maintain or increase the quality of the paper web 108 produced by the paper machine 102. In the following description, the machine direction (MD) of the web 108 denotes the direction along the (longer) length of the web 108.

In this example, the paper machine 102 includes at least one headbox 112, which distributes a pulp suspension uniformly across the machine onto a continuous moving wire screen or mesh 113. The pulp suspension entering the headbox 112 may contain, for example, 0.2-3% wood fibers, fillers, and/or other materials, with the remainder of the suspension being water.

Arrays of drainage elements 114, such as vacuum boxes, remove as much water as possible to initiate the formation of the web 108. An array of steam actuators 116 produces hot steam that penetrates the paper web 108 and releases the latent heat of the steam into the paper web 108. An array of rewet shower actuators 118 adds small droplets of water (which may be air atomized) onto the surface of the paper web 108. The paper web 108 is then often passed through a calender having several nips of counter-rotating rolls. Arrays of induction heating actuators 120 heat the shell surfaces of various ones of these rolls.

Two additional actuators 122-124 are shown in FIG. 1. A thick stock flow actuator 122 controls the consistency of incoming stock received at the headbox 112. A steam flow actuator 124 controls the amount of heat transferred to the paper web 108 from drying cylinders. The actuators 122-124 could, for example, represent valves controlling the flow of stock and steam, respectively. These actuators may be used for controlling the dry weight and moisture of the paper web 108. Additional flow actuators may be used to control the proportions of different types of pulp and filler material in the thick stock and to control the amounts of various additives (such as retention aid or dyes) that are mixed into the stock.

This represents a brief description of one type of paper machine 102 that may be used to produce a paper product. Additional details regarding this type of paper machine 102 are well-known in the art and are not needed for an understanding of this disclosure.

In order to control the paper-making process, one or more properties of the paper web 108 may be continuously or repeatedly measured. The web properties can be measured at one or various stages in the manufacturing process. This information may then be used to adjust the paper machine 102, such as by adjusting various actuators within the paper machine 102. This may help to compensate for any variations of the web properties from desired targets, which may help to ensure the quality of the web 108.

As shown in FIG. 1, the paper machine 102 includes one or more scanners 126-128, each of which may include one or more sensors. Each scanner 126-128 is capable of measuring one or more characteristics of the paper web 108. For example, each scanner 126-128 could include sensors for measuring the tension, caliper, moisture, anisotropy, basis weight, color, gloss, sheen, haze, surface features (such as roughness, topography, or orientation distributions of surface features), or any other or additional characteristics of the paper web 108.

Each scanner 126-128 includes any suitable structure or structures for measuring or detecting one or more characteristics of the paper web 108, such as one or more sets of sensors. The use of scanners represents one particular embodiment for measuring web properties. Other embodiments could be used, such as those including one or more stationary sets or arrays of sensors, deployed in one or a few locations across the web or deployed in a plurality of locations across the whole width of the web such that substantially the entire web width is measured.

The controller 104 receives measurement data from the scanners 126-128 and uses the data to control the paper machine 102. For example, the controller 104 may use the measurement data to adjust any of the actuators or other components of the paper machine 102. The controller 104 includes any suitable structure for controlling the operation of at least part of the paper machine 102, such as a computing device. Note that while a single controller 104 is shown here, multiple controllers 104 could be used, such as different controllers that control different variables of the web.

The network 106 is coupled to the controller 104 and various components of the paper machine 102 (such as the actuators and scanners). The network 106 facilitates communication between components of the system 100. The network 106 represents any suitable network or combination of networks facilitating communication between components in the system 100. The network 106 could, for example, represent a wired or wireless Ethernet network, an electrical signal network (such as a HART or FOUNDATION FIELDBUS network), a pneumatic control signal network, or any other or additional network(s).

The controller(s) 104 can operate to control one or more aspects of the paper machine 102 using one or more models. For example, each model could associate one manipulated variable with one controlled variable. A controlled variable generally represents a variable that can be measured or inferred and that is ideally controlled to be at or near a desired setpoint or within a desired range of values. A manipulated variable generally represents a variable that can be adjusted in order to alter one or more controlled variables.

In order to tune a controller 104, at least one operator console 130 can communicate with the controller 104 over a network 132. The operator console 130 generally represents a computing device that supports one or more techniques for robust tuning of MPC and other model-based controllers. The techniques for robust tuning of model-based controllers are described in more detail below. The network 132 represents any suitable network or combination of networks that can transport information, such as an Ethernet network.

In this example, the operator console 130 includes one or more processing devices 134, one or more memories 136, and one or more interfaces 138. Each processing device 134 includes any suitable processing or computing device, such as a microprocessor, microcontroller, digital signal processor, field programmable gate array, application specific integrated circuit, or discrete logic devices. Each memory 136 includes any suitable storage and retrieval device, such as a random access memory (RAM) or Flash or other read-only memory (ROM). Each interface 138 includes any suitable structure facilitating communication over a connection or network, such as a wired interface (like an Ethernet interface) or a wireless interface (like a radio frequency transceiver).

Note that while the operator console 130 is described as implementing the technique(s) for robust tuning of model-based controllers, other types of devices could also be used. For instance, the operator console 130 could interact with a server 140, and the server 140 could actually execute the algorithms used to implement one or more techniques for robust tuning of model-based controllers. In this case, the operator console 130 could present a graphical user interface and interact with a user. The server 140 could include one or more processing devices, one or more memories, and one or more interfaces (similar to the operator console 130).

Although FIG. 1 illustrates one example of a web manufacturing or processing system 100, various changes may be made to FIG. 1. For example, other systems could be used to produce other paper or non-paper products. Also, while shown as including a single paper machine 102 with various components and a single controller 104, the system 100 could include any number of paper machines or other machinery having any suitable structure, and the system 100 could include any number of controllers. In addition, FIG. 1 illustrates one example operational environment in which MPC or other model-based controller(s) can be tuned. This functionality could be used in any other suitable system.

In the following description, robust tuning techniques are described with respect to an MD-MPC tuning problem under the internal model control structure, which helps to simplify the tuning problem. “MD-MPC” here indicates that the controller being tuned is used to control an MD property of a web 108 using an MPC controller. Instead of tuning the weighting matrices in the MPC optimization problem, two filter parameters (which are referred to as λ-parameters) are used to adjust the closed-loop performance of the controller. It is assumed that a nominal process model is known, and parametric uncertainties on the process parameters are considered. Performance requirements are specified in terms of overshoot and settling time. This helps to maintain the friendliness of the proposed results to paper machine operators or end users but increases the difficulty of the analysis and parameter design. Due to the inevitable existence of time delays, analytical expressions for the closed-loop responses and the performance indices are not used. The structure of the resultant underlying optimization problem thus becomes unclear. Moreover, to obtain a satisfactory user experience, the computation time for the tuning procedure can be very limited.

Considering these difficulties, an efficient heuristic approach can be utilized to find a solution to the tuning problem. Example features of one heuristic approach include the following:

-   -   (1) Based on the small gain theorem, a sufficient condition for         robust stability under the parametric uncertainties is         presented.

(2) An efficient performance visualization technique is proposed, which allows the characterization of all possible step responses of a set of systems described by the parametric uncertainties. To improve the user experience, a fast implementation of the technique is also discussed.

(3) The parameter tuning problem is cast into a constrained optimization problem. To solve this problem, the empirical monotonic/unimodality properties of the overshoot and settling time with respect to the i-parameters are analyzed.

(4) Utilizing the visualization technique and the above properties, an iterative tuning algorithm is proposed to solve the optimization problem, based on which the model parameters can be tuned with a satisfactory performance within an acceptable computation time. The efficiency of the algorithm is demonstrated by an identified process model that is used for MD-MPC of a paper machine at an industrial site.

System Description and Problem Formulation

FIG. 2 illustrates an example internal model control structure 200 employed for MD-MPC according to this disclosure. This closed-loop system includes three main parts: an MD model 202, a MPC controller 204, and user-specified filters (F_(r) and F_(d)) 206-208.

In some embodiments, the MD model 202 represents a single-input, single-output (SISO) MD process 210 as a first-order-plus-dead-time (FOPDT) model, which can have the form:

$\begin{matrix} {{G_{0}(s)} = {\frac{g_{0}}{{T_{p\; 0}s} + 1}^{{- t_{d\; 0}}s}}} & (1) \end{matrix}$

where g₀, T_(p0), and t_(d0) indicate the process gain, time constant, and time delay, respectively, of the MD process. A discrete-time state-space realization of the model in Equation (1) can be written as:

$\begin{matrix} \left\{ {\begin{matrix} {{{x\left( {k + 1} \right)} = {{{Ax}(k)} + {{Bu}(k)}}},} \\ {{y(k)} = {{{Cx}(k)} + {d(k)}}} \end{matrix}.} \right. & (2) \end{matrix}$

Note that the time delay of the model in Equation (1) has been absorbed in (A, B, C, D) so that u(k) can be used instead of u(k−T_(d0)) (T_(d0) is a discretized version of t_(d0)). Since using Δu(k) as an input simplifies the predictions in MPC iterations, the state-space model in Equation (2) can be restructured as follows:

$\begin{matrix} {{\begin{bmatrix} {\Delta \; {x\left( {k + 1} \right)}} \\ {{Cx}\left( {k + 1} \right)} \end{bmatrix} = {{\underset{\underset{A_{a}}{}}{\begin{bmatrix} A & 0 \\ {CA} & 1 \end{bmatrix}}\underset{\underset{x_{a}{(k)}}{}}{\begin{bmatrix} {\Delta \; {x(k)}} \\ {{Cx}(k)} \end{bmatrix}}} + {\underset{\underset{B_{a}}{}}{\begin{bmatrix} B \\ {CB} \end{bmatrix}}\Delta \; {u(k)}}}},{{y(k)} = {{\underset{\underset{C_{a}}{}}{\begin{bmatrix} 0 & 1 \end{bmatrix}}{x_{a}(k)}} + {{d(k)}.}}}} & (3) \end{matrix}$

which is used in the derivation of the MD-MPC solution.

To generate an MPC solution (considering the MPC for MD control), a prediction model can be used to obtain estimations of controlled variables in a time horizon. Let H_(u) denote the control horizon and H_(p) denote the prediction horizon (where 1≦H_(u)≦H_(p)). Based on Equation (3), the prediction model can be derived as follows:

$\begin{matrix} \left\{ {\begin{matrix} {{{\hat{x}}_{a}\left( {k + i} \right)} = \begin{matrix} {{A_{a}^{i}{{\hat{x}}_{a}(k)}} + \sum_{j = 1}^{m\; i\; n{\{{H_{u},i}\}}}} \\ {{A_{a}^{i - j}B_{a}\Delta \; {u\left( {k + j - 1} \right)}},} \end{matrix}} \\ {{\hat{y}\left( {k + i} \right)} = {C_{a\;}{{\hat{x}}_{a}\left( {k + i} \right)}}} \end{matrix},} \right. & (4) \end{matrix}$

where i=1, 2, . . . , H_(p). In MD-MPC, the following cost function can be defined for obtaining a control move:

$\begin{matrix} {{{\min\limits_{\Delta \; U}} = {{\left( {\hat{Y} - Y_{ref}} \right)^{T}{Q_{1}\left( {\hat{Y} - Y_{ref}} \right)}} + {\Delta \; U^{T}Q_{2}\Delta \; U} + {\left( {U - U_{ref}} \right)^{T}{Q_{3}\left( {U - U_{ref}} \right)}}}},\mspace{20mu} {{where}\text{:}}} & (5) \\ {\mspace{20mu} {{\hat{Y} = \begin{bmatrix} {\hat{y}\left( {k + 1} \right)} \\ {\hat{y}\left( {k + 2} \right)} \\ \vdots \\ {\hat{y}\left( {k + H_{p}} \right)} \end{bmatrix}},{{\Delta \; U} = \begin{bmatrix} {\Delta \; {u(k)}} \\ {\Delta \; {u\left( {k + 1} \right)}} \\ \vdots \\ {\Delta \; {u\left( {k + H_{u} - 1} \right)}} \end{bmatrix}},\mspace{20mu} {U = {{\begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{bmatrix}{u\left( {k - 1} \right)}} + {\begin{bmatrix} 1 & 0 & \ldots & 0 \\ 1 & 1 & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 1 & \ldots & 1 & 1 \end{bmatrix}\Delta \; U}}},}} & (6) \end{matrix}$

Here, U_(ref) and Y_(ref) are dimension-compatible vectors containing reference signals of the corresponding variables at different time instants. Also, Q₁, Q₂, and Q₃ are weighting matrices.

In practice, physical constraints on actuators can be taken into account when solving the above quadratic programming (QP) problem. However, it is well-known that the solution of a constrained QP problem may need to be numerically computed, and this process can take a long time (which can be undesirable in various circumstances). To simplify the tuning analysis, an unconstrained MPC problem can be considered instead. Thus, a closed-form solution can be easily derived:

$\begin{matrix} \begin{matrix} {{\Delta \; {u(k)}} = {f\left( {U_{ref},Y_{ref}} \right)}} \\ {= {{{z\left\lbrack {\left( {z - 1} \right) - {\left( {z - 1} \right)\xi \; K_{x}{\theta (z)}} - {\xi \; K_{u - 1}}} \right\rbrack}^{- 1}\xi \; K_{yt}Y_{ref}} +}} \\ {{{z\left\lbrack {\left( {z - 1} \right) - {\left( {z - 1} \right)\xi \; K_{x}{\theta (z)}} - {\xi \; K_{u - 1}}} \right\rbrack}^{- 1}\xi \; K_{xt}{U_{ref}.}}} \end{matrix} & (7) \end{matrix}$

The expressions of ξ, K_(x), K_(u-1), K_(xt), and K_(yt) are presented below and can be obtained by standard derivations in MPC.

With respect to the two user-specified filters F_(r) and F_(d) 206-208 in FIG. 2, these filters 206-208 can significantly affect the controller performance. These filters 206-208 are used respectively for filtering the output target y_(tgt)(k) and the estimated disturbance {circumflex over (d)}(k). Thus, a two-degrees of freedom (2DOF) model predictive controller can have two main tuning parameters, namely (i) the desired closed-loop time constant for disturbance rejection and (ii) the desired closed-loop time constant for setpoint tracking. With the filtered signals, a reference trajectory can be calculated as follows:

$\begin{matrix} {{Y_{ref}(k)} = {\begin{bmatrix} {y_{ref}\left( {k + 1} \right)} \\ \vdots \\ {y_{ref}\left( {k + H_{p}} \right)} \end{bmatrix} = {{F_{r}{y_{tgt}(k)}} = {F_{d}{{\hat{d}(k)}.}}}}} & (8) \end{matrix}$

Here, y_(ref)(k)=f_(r)(z)y_(tgt)(k)−f_(d)(z){circumflex over (d)}(k), where f_(r)(z) and f_(d)(z) are a reference tracking filter and a disturbance rejecting filter, which can be defined as follows:

$\begin{matrix} {{{f_{r}(z)} = {\frac{b_{r}z^{- 1}}{1 - {a_{r}z^{- 1}}}z^{- T_{d}}}},{{f_{d}(z)} = {\frac{b_{d}z^{- 1}}{1 - {a_{d}z^{- 1}}}z^{- T_{d}}}},} & (9) \end{matrix}$

where

${a_{r} = ^{- \frac{T_{s}}{\lambda \; T_{p}}}},$

b_(r)=1−a_(r),

${a_{d} = ^{- \frac{T_{s}}{\lambda_{d}T_{p}}}},$

b_(d)=1−a_(d). Parameters λ λ_(d) are called the reference tracking performance ratio and the disturbance rejecting performance ratio, respectively. From the state-spaces of f_(r)(z) and f_(d)(z), F_(r) and F_(d) can be readily constructed, such as by using a procedure similar to the construction of the prediction model in Equation (4).

Note that the two filters F_(r) and F_(d) 206-208 enable one to treat reference tracking and disturbance rejection separately. Moreover, through the design of the filter F_(r) 206, the response of the output to the target signal can be directly controlled without affecting the disturbance and noise rejection (if it exists). With the help of this 2DOF control, some user-familiar features of the system (such as overshoot, settling time, or tracking error) can be linked to λ and λ_(d), which simplifies the following closed-loop tuning and makes tuning analysis more user-friendly.

Additional details regarding the use of the F_(r) and F_(d) filters 206-208 to obtain 2DOF are provided in U.S. patent application Ser. No. 13/907,495 filed on May 31, 2013, which is hereby incorporated by reference in its entirety.

The MD-MPC tuning problem can therefore involve determining how to manipulate λ and λ_(d) so that (i) the closed-loop system in FIG. 2 is robustly stable and (ii) one or more controlled variables will track one or more output targets with small overshoot and quick response. Although there are advantages to using 2DOF control, it cannot be denied that tuning two parameters is more difficult than tuning one parameter. As a result, the control structure in FIG. 2 can be revisited, and a check can be made whether the tuning of λ and λ_(d) can be separated. Looking at the stability of the closed-loop system (as it could be of the highest priority during tuning), it can be seen from FIG. 2 and the expression of F_(r) that the system's stability is only dependent on λ_(d). Therefore, it is possible to ignore the effects of λ when examining the closed-loop stability. Moreover, it can be seen that the main effect of λ on F_(r) is to change the time constant of the filter F_(r) 206, which plays a role in balancing overshoot with response speed. In other words, as λ increases, the overshoot of the closed-loop system becomes smaller, while the system's response gets slower.

Robust Stability Analysis

This section discusses the tuning of λ and Δ_(d). Robust stability analysis, as one part of tuning λ_(d), is discussed first. The outcome of the analysis in this section results in a feasible region λ_(d)ε[λ*_(d),∞] such that the closed-loop system is stable, which provides reference for an iterative parameter tuning procedure described later.

FIGS. 3 through 12 illustrate details of an example technique for specifying and visualizing robust tuning of model-based controllers according to this disclosure. For simplicity, the closed-loop system of FIG. 2 is rearranged as shown in FIG. 3. In FIG. 3, the input U_(ref) and d have been removed as they do not change the stability of the whole system. Moreover, the MD model 202 has been combined with the F_(d) filter 208 for brevity. This is indicated as F_(d) =[F_(d) ¹ ,F_(d) ²], where F_(d) ¹ and F_(d) ² are the transfer functions from y to Y_(ref) and from y to Y_(ref); respectively. In FIG. 3, the actual MD process 210 is represented by P, which is different from the MD model 202 and can be written in the multiplicative uncertainty form as follows:

P=G ₀(1+WΔ)  (10)

where Δ indicates the model uncertainty and W is a weighting function.

By pulling out Δ, the closed-loop system can be represented in the general structure 400 (see FIG. 4A). The expression of N and the transfer function of the closed-loop system in FIG. 4A can be expressed as follows:

$\begin{matrix} {N = \begin{bmatrix} N_{11} & N_{12} \\ N_{21} & N_{22} \end{bmatrix}} & (11) \\ {{F\left( {N,\Delta} \right)} = {N_{22} + {N_{21}{\Delta \left( {I - {N_{11}\Delta}} \right)}^{- 1}{N_{12}.}}}} & (12) \end{matrix}$

It is known that the stability of the closed-loop system F(N,Δ) in FIG. 4A is equivalent to the stability of the system 402 in FIG. 4B (M=N₁₁). Therefore, the M−Δ system is examined instead of the whole system. In the MD-MPC system:

M=−(1+F ₂ F _(d) ¹ +F ₂ F _(d) ² G ₀)⁻¹ F ₂ F _(d) ² G ₀ W.  (13)

For the M−Δ system, the robust stability has been well investigated. This system may be robustly stable if and only if:

det(1−M(jω)Δ)≠0,∀ω,∀|Δ|≦1.  (14)

For SISO systems, the condition in Equation (14) is equivalent to:

|M(jω)<1,∀ω,∀|Δ|≦1

|W(jω)T(jω|<1,∀w,  (15)

where

$T = {\frac{b_{d}z^{- 1}}{1 - {a_{d}z^{- 1}}}z^{- T_{d}}}$

is the sensitivity function of the MD-MPC system.

Using the condition in Equation (15), a feasible region of λ_(d) can be obtained such that robust stability is guaranteed. In this procedure, the weighting function W can be chosen so that the obtained region of λ_(d) is not too conservative. From the expression of P, Δ=W⁻¹(P−G₀)G₀ ⁻¹. As a part of the robust stability condition, |Δ|≦1, which implies that:

|W ⁻¹(P−G ₀)G ₀ ⁻¹|≦1 and |W|≧|(P−G ₀)G ₀−|.  (16)

This provides a way to construct the weighting function W using the maximum value of |(P−G₀)G₀ ⁻¹|, which is known as the multiplicative error.

It is assumed here that the model uncertainty is only reflected in the model parameters because parametric uncertainty is relatively easy to understand for users, even for users without control backgrounds. With this consideration, for a FOPDT model with the following parametric uncertainty:

g ₀ −Δg≦g≦g ₀ +Δg,

T _(p0) −ΔT _(p) ≦T _(p) ≦T _(p0) ≦ΔT _(p),

t _(d0) −Δt _(d) ≦t _(d) ≦t _(d0) ×Δt _(d),  (17)

the multiplicative error can be expressed as:

$\begin{matrix} {{\left( {P - G_{0}} \right)G_{0}^{- 1}} = {\frac{{g\; {^{{({t_{d\; 0} - t_{d}})}s}\left( {{T_{p\; 0}s} + 1} \right)}} - {g_{0}\left( {{T_{p}s} + 1} \right)}}{g_{0}\left( {{T_{p}s} + 1} \right)}.}} & (18) \end{matrix}$

The weighting function W can be constructed in any suitable manner, such as by following a similar argument as that in Hu et al., “Systematic h_(∞) weighting function selection and its application to the real-time control of a vertical take-off aircraft,” Control Engineering Practice, vol. 8, pp. 241-252, 2000 (which is hereby incorporated by reference in its entirety). The construction details and the expression of W(s) are shown below. Note that W(s) plays the role of translating parametric uncertainty into multiplicative uncertainty, which facilitates the tuning of λ_(d) with classic robust control theory.

To determine the feasible region of λ_(d), one approach is to find all λ_(d)'s such that |T(jω)| is less than |W(jω)| over all frequencies. From the expression of T(s), the following can be obtained:

$\begin{matrix} {{T(s)} = {{\frac{1}{{\lambda_{d}T_{p\; 0}s} + 1}^{{- t_{d\; 0}}s}\mspace{14mu} {and}\mspace{14mu} {{T(\omega)}}} = {\frac{1}{\sqrt{\left( {\lambda_{d}T_{p\; 0}\omega} \right)^{2} + 1}}.}}} & (19) \end{matrix}$

Thus, increasing λ_(d) reduces |T(jω)| as well as the bandwidth of T(s) (see FIG. 5 containing a graph 500 showing the effects of different λ_(d) values). This means that the system becomes more robust as λ_(d) increases. As a result, to get the feasible region of λ_(d), the minimum λ_(d) (denoted λ*_(d)) can be identified so that the robust stability condition is satisfied.

Performance Visualization & Fast Implementation

In this section, the performance visualization problem used in the parameter tuning procedure is considered. One objective of the visualization can be to graphically characterize the envelopes of the responses of a set of systems whose parameters lie within user-specified uncertainty intervals, subject to a step setpoint or load disturbance change, given the values of λ and λ_(d). This can be useful in industrial applications as it allows operators to easily tell the resultant consequences of choosing a combination of λ and λ_(d).

Mathematically, the performance visualization problem can be viewed as the composition of two sequences of optimization problems:

Problem 1: For all t=1, 2, . . . , solve respectively

$\max\limits_{g,t_{d},T_{p}}{y(t)}$ ${{s.t.\mspace{14mu} t_{d}} \in \left\lbrack {{\underset{\_}{t}}_{d},{\overset{\_}{t}}_{d}} \right\rbrack};$ ${T_{p} \in \left\lbrack {{\underset{\_}{T}}_{p},{\overset{\_}{T}}_{p}} \right\rbrack};$ ${g \in \left\lbrack {\underset{\_}{g},\overset{\_}{g}} \right\rbrack};$ X_(a)(k + 1) = A_(a)X_(a)(k) + B_(a)Δ u(k), y(k) = C_(a)X_(k) + d(k), Δ u(k) = f_(MPC)(U_(tgt), Y_(tgt), X_(a)(k)), k = 1, 2, …  , t; X_(a)(0) = 0; and $\min\limits_{g,t_{d},T_{p}}{y(t)}$ ${{s.t.\mspace{14mu} t_{d}} \in \left\lbrack {{\underset{\_}{t}}_{d},{\overset{\_}{t}}_{d}} \right\rbrack};$ ${T_{p} \in \left\lbrack {{\underset{\_}{T}}_{p},{\overset{\_}{T}}_{p}} \right\rbrack};$ ${g \in \left\lbrack {\underset{\_}{g},\overset{\_}{g}} \right\rbrack};$ X_(a)(k + 1) = A_(a)X_(a)(k) + B_(a)Δ u(k), y(k) = C_(a)X_(a)(k) + d(k), Δ u(k) = f_(MPC)(U_(tgt), Y_(tgt), X_(a)(k)), k = 1, 2, …  , t; X_(a)(0) = 0,

where f_(MPC)(U_(tgt), Y_(tgt), X_(a)(k)) denotes the optimal solution to the constrained MPC problem.

Note that model mismatch is considered in the above optimization problems since Δu(k) is calculated according to the nominal process parameters. Several issues exist in finding the exact solutions to Problem 1, including:

(1) y(t) is not necessarily a convex function of the optimization parameters g, T_(p), and t_(d) (in fact, g and T_(p) are explicitly expressed in A_(a), while t_(d) affects A_(a) implicitly by controlling the size of the matrix);

(2) the complexity of the dependence of y(t) on g, T_(p), and t_(d) increases with an increase of t; and

(3) as performance visualization is one step in the overall iterative tuning procedure, the computation time allowed to solve the above problem can be very limited.

In light of these issues, efficient low-complexity heuristic solutions to this problem can be used. For example, Problem 1 could be rewritten into an equivalent form:

Problem 2: For all t=1, 2, . . . , solve respectively

$\max\limits_{g,t_{d},T_{p}}{h\left( {t,g,t_{d},T_{p},U_{tgt},Y_{tgt},{X_{a}(0)},d} \right)}$ ${{s.t.\mspace{14mu} t_{d}} \in \left\lbrack {{\underset{\_}{t}}_{d},{\overset{\_}{t}}_{d}} \right\rbrack},{T_{p} \in \left\lbrack {{\underset{\_}{T}}_{p},{\overset{\_}{T}}_{p}} \right\rbrack}$ ${g \in \left\lbrack {\underset{\_}{g},\overset{\_}{g}} \right\rbrack},{and}$ $\min\limits_{g,t_{d},T_{p}}{h\left( {t,g,t_{d},T_{p},U_{tgt},Y_{tgt},{X_{a}(0)},d} \right)}$ ${{s.t.\mspace{14mu} t_{d}} \in \left\lbrack {{\underset{\_}{t}}_{d},{\overset{\_}{t}}_{d}} \right\rbrack},{T_{p} \in {p\left\lbrack {{\underset{\_}{T}}_{p},{\overset{\_}{T}}_{p}} \right\rbrack}},{g \in \left\lbrack {\underset{\_}{g},\overset{\_}{g}} \right\rbrack},$

where h(t, g, t_(d), T_(p), U_(tgt), Y_(tgt), X_(a)(0), d) is obtained by composing y(t) with the state-space equation and optimal MPC law f_(MPC)(U_(tgt), Y_(tgt), X_(a)(k)).

In Problem 2, it can be observed that both optimization problems are nonlinear optimization problems within polytopes. According to the Karush-Kuhn-Tucker (KKT) necessary conditions, an optimizer can be obtained by enumerating all possible combinations of active constraints and solving the resultant unconstrained problems. As a result, one efficient heuristic can be to assume that the optimizer is achieved at one of the vertices of the polytopes, and the problems can therefore be solved by comparing the values of the objective function for only 2³ points in the parameter space, which results in the following algorithm.

Algorithm 1 Performance visualization procedure 1: Calculate the required responses y_(i) for the vertex i of the parametric set (i = 1, 2, . . . , 8) for time instants t = 1, 2, . . . , N; 2: Obtain the upper envelope y(t) by calculating max_(i∈{1, 2, . . . , 8})y_(i)(t) for t = 1, 2, . . . , N; 3: Obtain the lower envelope y(t) by calculating min_(i∈{1, 2, . . . , 8})y_(i)(t) for t = 1, 2, . . . , N. Although this heuristic is not guaranteed to be optimal, it is intuitive in that extreme behaviors of the step responses mostly happen at the extreme process parameters. Also, this method is applicable to the characterization of envelopes for other process variables, such as control variables (see below for details).

One potential benefit of the above-proposed performance visualization technique is that, for a user-specified set of processes, it allows direct evaluation of time-domain performance indices, the expressions of which are normally not possible to calculate analytically. It also provides feedback information for the overall iterative parameter tuning procedure.

Assuming all responses in a set have the same final values (and this assumption is automatically satisfied since the system under consideration is of type 1), the definitions of two well-accepted performance indices for the set of step responses are introduced, which are generalized from their counterparts for a single response.

Definition 1 (Overshoot): The overshoot OS of a set of step responses with the same final value is the maximum value in all responses minus the final value divided by the final value.

Definition 2 (Settling time): The settling time t_(s) of a set of step responses with the same final value is the time required for all responses to reach and stay within a range of pre-specified percentage of the final value (a ±5% value is assumed below, although other values could be used).

Based on these definitions and the proposed visualization method, the values of OS and t_(s) can be calculated numerically. For example, the OS can be computed by finding the maximum peak in y(t) and:

$\begin{matrix} {{OS} = \frac{{\max \left\{ {\overset{\_}{y}(t)} \right\}} - {y(\infty)}}{y(\infty)}} & (20) \end{matrix}$

The settling time can be calculated by reversing the time series y and y and identifying the time tag of the first element that escapes the ±5% interval.

One way to implement the proposed procedure is to run a real-time simulation using the HONEYWELL MPC and simulator, which is relatively time consuming (compared with the tuning procedure) and adds to the computational burden. In this regard, a fast simulator that only considers the unconstrained MPC (which reduces to an LQR controller) can be embedded in the visualization procedure. The fast simulator can be implemented in state-space models, which are compatible with both SISO and multiple-input, multiple-output (MIMO) cases. Satisfactory performance can be achieved for most cases (see below for details).

The visualization technique can be illustrated using extensive simulations. In the following, the proposed procedure is applied to FOPDT processes, and the results are shown in FIGS. 6A-6I. Three types of typical processes are considered. The first one has balanced time constant and delay (FIGS. 6A-6C), the second one is delay dominant (FIGS. 6D-6F), and the third one is time-constant dominant (FIGS. 6G-6I). For each type of process, the values of (λ, λ_(d)) are set to (1.5,1.5), (2,1), and (1,2), respectively. To demonstrate the efficiency of the visualization procedure, the responses of one hundred randomly generated systems satisfying the parameter uncertainty levels are also plotted in each case (the uncertainty levels are taken as ±24% of the nominal process parameters). The parameters of each case and the measured overshoots and settling times are shown in FIGS. 6A-6I.

From FIGS. 6A-6I, it is observed that the resultant envelopes 602 a-602 i may efficiently characterize sets of possible responses generated by the processes with parameters lying within the specified intervals, even when the systems are unstable (such as in FIG. 6E). The numerical evaluations of settling times and overshoots are also shown to be accurate. Notice that the computation complexity of the procedure is satisfactory, as a single run of the procedure for a FOPDT process takes only about 0.17 seconds on a laptop with an INTEL CORE-I5 processor and 6 GB of memory.

Iterative Parameter Tuning

In this section, an iterative tuning procedure is proposed for λ and λ_(d) based on the stabilizing region of λ determined in the “ROBUST STABILITY ANALYSIS” section and the proposed visualization techniques in the “PERFORMANCE VISUALIZATION & FAST IMPLEMENTATION” section. For notational simplicity, t_(s)(λ, λ_(d)) and OS(λ, λ_(d)) are used to represent the relationship of the settling time and the overshoot with λ and λ_(d).

Before proceeding further, example performance requirements for parameter tuning are first presented, and the tuning problem is formulated into a constrained optimization problem. Here, time-domain specifications are employed to quantify the performance of a system, such as overshoot and settling time. In industrial applications, one possible preference of users on overshoot is that overshoot should not exceed a certain level. Since overshoot is a unified variable that does not depend on the time constant or time delay of the system, this preference can be easily implemented as a constraint that restricts the overshoot within a specified region. A smaller settling time may be more preferable, but the value of the settling time implicitly depends on the time constant and time delay. Therefore, in some embodiments, the settling time could be minimized in the tuning procedure. In this way, the parameter tuning problem can be formulated into the following optimization problem:

$\begin{matrix} {{\min\limits_{\lambda,\lambda_{d}}{t_{s}\left( {\lambda,\lambda_{d}} \right)}}{{s.t.\mspace{14mu} {{OS}\left( {\lambda,\lambda_{d}} \right)}} \leq {{OS}^{*}.}}} & (21) \end{matrix}$

In order to develop algorithms to find pairs (λ, λ_(d)) that solve Equation (21), geometric properties of t_(s)(λ, λ_(d)) and OS(λ, λ_(d)) can be explored. However, analytical expressions of t_(s)(λ, λ_(d)) and OS(λ, λ_(d)) in general do not exist, especially for discrete-time delayed cases. With this in mind, one approach that can be taken is to provide a qualitative analysis on the properties of the function and then to numerically verify these properties.

As shown above in the “ROBUST STABILITY ANALYSIS” section, the stability of the system is determined by λ_(d), and the robustness of the system increases with increases of λ_(d). It is therefore intuitive from an engineering perspective that t_(s)(λ, λ_(d)) is a unimodal function of λ_(d) for a fixed value of A. A very small value of λ_(d) causes a large settling time due to relative aggressive and oscillatory responses, while a large value of λ_(d) also causes a large settling time due to over-sluggish responses. Similarly, when λ_(d) is fixed, t_(s)(λ, λ_(d)) is a unimodal function of A. On the other hand, the filter F_(r) 206 controls only the speed of the response, and thus a larger value of λ leads to a smaller value of overshoot. In this way, OS(λ, λ_(d)) can be empirically treated as a monotonically decreasing function of λ.

To verify this analysis, the relationship of t_(s)(λ, λ_(d)) and OS(λ, λ_(d)) with λ and λ_(d) for an FOPDT process can be numerically evaluated (see FIGS. 7A and 7B where plots 700-702 respectively display overshoot and settling time against λ and λ_(d) values). The envelopes in the visualization procedure can be generated within T_(p)ε[(1−r)T_(p0),(1+r)T_(p0)], t_(d)ε[(1−r)t_(d0),(1+r)t_(d0)], gε[(1−r)_(g0),(1+r)g₀], r being the uncertainty level. To improve visibility, large values of overshoot and settling time are truncated without affecting the verification results. Instead of the value of overshoot, the logarithm of overshoot is plotted to verify the monotonicity property of OS(λ, λ_(d)). From FIGS. 7A and 7B, the monotonic property of OS(λ, λ_(d)) with respect to λ and the unimodal property of t_(s)(λ, λ_(d)) with respect to λ_(d) are observed, which verifies the empirical analysis.

Based on the empirical properties of t_(s)(λ, λ_(d)) and OS(λ, λ_(d)), an efficient and robust tuning algorithm is proposed below. The efficient and robust properties of the algorithm are emphasized since:

(i) In practice, a satisfactory solution obtained within a small amount of time results in a better user experience, compared to that of an optimal solution obtained based on certain models (with unavoidable modeling errors) at the cost of more computation time; and

(ii) The explicit expressions of t_(s)(λ, λ_(d)) and OS(λ, λ_(d)) are unknown and only numerical values are available using the visualization techniques, which limits the allowed types of algorithms that can be considered. For instance, Newton and quasi-Newton algorithms may be prohibited since the numerical evaluation of the first-order derivatives could lead to unexpected errors and thus the failure of the overall tuning procedure.

In light of these factors, in some embodiments, line-search methods are used to find a satisfactory combination of λ and λ_(d). In this approach, the algorithm is performed iteratively to find a pair of λ and λ_(d) values that heuristically solves Equation (21).

Algorithm 2 Tuning of λ and λ_(d)  1: Input the uncertainty intervals [T _(p), T _(p)], [t _(d), t _(d)] and [g, g] for process parameters T_(p0), t_(d0) and g₀, respectively;  2: Input the overshoot specification OS*;  3: Input ε and N;  4: λ _(d) ← λ _(d)*; λ _(d) ← 100;  5: λ ← 0.02; λ_(d) ← λ _(d); i ← 0;  6: λ* ← λ; λ_(d)* ← λ _(d); t_(s)* = ( T _(p) + t _(d)) × 10; t_(s) ⁰ ← 0;  7: while t_(s) ⁰ ≠ t_(s)(λ, λ_(d)) && i < N do  8: t_(s) ⁰ = t_(s)(λ, λ_(d));  9: i ← i + 1;

 Stage 1: tuning of λ_(d) 10: λ_(d1) ← λ _(d) + ( λ _(d) − λ _(d)) × 0.382; 11: λ_(d2) ← λ _(d) + ( λ _(d) − λ _(d)) × 0.618; 12: while λ_(d) − λ _(d) > ε do 13: if t_(s)(λ,λ_(d1)) > t_(s)(λ,λ_(d2)) then 14: λ _(d) ← λ_(d1); λ_(d1) ← λ_(d2); 15: λ_(d2) ← λ _(d) + ( λ _(d) − λ _(d)) × 0.618; 16: else 17: λ _(d) ← λ_(d2); λ_(d2) ← λ_(d1); 18: λ_(d1) ← λ _(d) + ( λ _(d) − λ _(d)) × 0.382; 19: end if 20: end while 21: λ_(d) ← ( λ _(d) + λ _(d))/2; 22: λ ← 0.02; λ ← 100;

 Stage 2: finding a proper λ 23: while λ − λ > ε do 24: λ ← λ + ( λ − λ) × 0.5; 25: if OS(λ,λ_(d)) − OS* > 0 then 26: λ ← λ; 27: else 28: λ ← λ; 29: end if 30: end while 31: λ ← 100;

 Stage 3: tuning of λ 32: λ₁ ← λ + ( λ − λ) × 0.382; 33: λ₂ ← λ + ( λ − λ) × 0.618; 34: while λ − λ > ε do 35: if t_(s)(λ₁,λ_(d)) > t_(s)(λ₂,λ_(d)) then 36: λ ← λ₁; λ₁ ← λ₂; 37: λ₂ ← λ + ( λ − λ) × 0.618; 38: else 39: λ ← λ₂; λ₂ ← λ₁; 40: λ₁ ← λ + ( λ − λ) × 0.382; 41: end if 42: end while 43: λ ← ( λ + λ)/2; 44: if t_(s)(λ,λ_(d)) < t_(s)* then 45: t_(s)* ← t_(s)(λ,λ_(d)); λ* ← λ; λ_(d)* ← λ_(d); 46: end if 47: end while 48: Output λ* and λ_(d)*; 49: end

Each iteration of this algorithm includes three stages, which are described in more detail below. The algorithm stops when a stationary point is achieved (meaning the same pair of λ and λ_(d) values is obtained in consecutive iterations) or when the computation time runs out (which could be counted as the number of iterations allowed). The algorithm collects the feasible (λ, λ_(d)) pair with the smallest settling time during the iterative procedure. This assumes that the nominal process parameters T_(p0), t_(d0), and g₀ are known, which can be readily obtained (such as by using the HONEYWELL PROFIT DESIGN STUDIO). Note that to avoid the possible conservativeness of the small gain theorem, the initial value of λ _(d) can be chosen as a small value (say 0.02), instead of λ*_(d) (see line 4).

In the first stage of the algorithm (lines 10-21), λ_(d) is first tuned to obtain the optimal settling time for a fixed λ. At this stage, λ is set to a very small value, which would lead to a relatively aggressive response for a fixed λ_(d). This choice simplifies the tuning procedure in that, when the tuning of λ_(d) is completed, it suffices to increase λ to find λ, that activates the overshoot constraint (stage 2) by the empirical monotonicity. By the empirical unimodality property, the optimization of the settling time with respect to λ_(d) is achieved by a line-search with almost a linear convergence rate without requiring numerically calculating the derivatives (a guaranteed linear convergence rate can be achieved, such as by using Fibonacci sequences).

In the second stage (lines 17-29), λ is further tuned to find λ that activates the constraint on overshoot (namely, the constraint becomes an equality constraint), based on the tuned value of λ_(d) in Stage 1. To do this, the following optimization problem can be considered:

$\begin{matrix} {{\min\limits_{\lambda}\left\lbrack {{{OS}\left( {\lambda,\lambda_{d}} \right)} - {OS}^{*}} \right\rbrack^{2}}{{{s.t.\mspace{14mu} t_{d}} \in \left\lbrack {{\underset{\_}{t}}_{d},{\overset{\_}{t}}_{d}} \right\rbrack},{T_{p} \in \left\lbrack {{\underset{\_}{T}}_{p},{\overset{\_}{T}}_{p}} \right\rbrack},{g \in {\left\lbrack {\underset{\_}{g},\overset{\_}{g}} \right\rbrack.}}}} & (22) \end{matrix}$

From the empirical monotonicity property of the overshoot with respect to λ, the objective function [OS(λ,λ_(d))−OS*]² is a unimodal function with respect to λ. Thus, Equation (22) can be solved by line-search algorithms, which leads to the codes in lines 22-30.

In the third stage (lines 31-42), the settling time is further optimized with respect to λ within the region λ≧λ, where λ is calculated in the second stage and λ_(d) is calculated in the first stage. Notice that, due to the monotonicity property of OS(λ, λ_(d)) with respect to λ, the overshoot constraint in Equation (21) is satisfied for all λ≧λ. Within this region, OS(λ, λ_(d)) is either a unimodal function or a monotonic function of A, which allows line-search algorithms to find the λ that achieves the smallest settling time.

Note that although optimization procedures are iteratively employed in Algorithm 2, the tuning result may not be optimal in either overshoot or settling time. As mentioned before, what is actually achieved is an efficient, fast, and robust tuning procedure that yields a combination of λ and λ_(d) that approximately solves Equation (21), resulting in satisfactory performance to end users (which can be jointly guaranteed by Algorithms 1 and 2).

Industrial Example

In this section, the efficiency of the proposed tuning algorithm is illustrated with the tuning results obtained from a typical SISO process in a paper machine. Consider the following model for MD-MPC control:

$\begin{matrix} {{G_{0}(s)} = {\frac{0.0135}{{60\; s} + 1}^{{- 90}\; s}}} & (23) \end{matrix}$

For this model, the proposed algorithm is applied for different levels of uncertainty and different specifications on overshoot. To test the performance of the algorithm in terms of the optimality of the settling time, a brutal search is performed for each level of uncertainty and overshoot specification.

The performance comparison is presented in FIGS. 8 through 10, where the computational time of the proposed algorithm is also indicated for each point. In FIG. 8, a graph 800 plots the overshoot and settling time obtained using the brutal search (line 802) and using the described algorithm (line 804) for one set of parametric uncertainties on the process gain, time constant, and time delay. In FIG. 9, a graph 900 plots the overshoot and settling time obtained using the brutal search (line 902) and using the described algorithm (line 904) for another set of parametric uncertainties on the process gain, time constant, and time delay. In FIG. 10, a graph 1000 plots the overshoot and settling time obtained using the brutal search (line 1002) and using the described algorithm (line 1004) for a third set of parametric uncertainties on the process gain, time constant, and time delay.

Although the outcome of the proposed algorithm does not have guaranteed optimality, it is consistently close to the result of the brutal search. The brutal search could ordinarily take around ten minutes to calculate a (λ, λ_(d)) pair for a given specification of overshoot by enumerating all points over a pre-specified gridded parameter region. On the other hand, the computation time of the proposed algorithms is satisfactory, taking around seven to fifteen seconds.

To take a closer look at the tuning parameters and the performance indices, the tuning results for uncertainty level—50%˜100% are presented in Table I below, and the closed-loop step response for uncertainty level—50%˜100% and OS*=20% is shown in FIG. 11.

TABLE I Tuning results for uncertainty level −50%~100% Overshoot specification OS* 10% 20% 30% 40% 50% Proposed OS    10% 16.98% 29.90% 33.42% 44.6159% algorithm t_(s) 2505 s 2415 s 2340 s 2340 s 2355 s λ 8.8279 8.0097 6.6243 6.2889 5.3149 λ_(d) 3.61  3.7413 3.9537 4.0039 4.1662 Brutal OS 9.1214% 16.98% 29.90% 30.3078% 44.6159% search t_(s) 2535 s 2415 s 2340 s 2340 s 2355 s λ 9.0214 8.0097 6.6243 6.5981 5.3149 λ_(d) 3.62  3.7413 3.9537 4.095 4.1662 In FIG. 11, a line 1102 represents the closed-loop step response of a system, and lines 1104-1106 define an envelope around the closed-loop step response (where the envelope is based on the parametric uncertainties).

To test the design results in a real-time MPC plus simulator environment, to account for model mismatch, the actual process is taken as:

$\begin{matrix} {{G(s)} = {\frac{0.0108}{{90\; s} + 1}^{{- 120}\; s}}} & (24) \end{matrix}$

which lies within the uncertainty region. The system is discretized at sampling time T_(s)=10 s. The MPC weighting matrices are set to Q₁=1, Q₂=0, and Q₃=0, the prediction horizon H_(p) is 42, and the control horizon H_(u) is 20. The constraints are taken as:

1895≦U≦6064

−379≦ΔU≦379  (25)

The initial operating point is y(0)=432 and u(0)=3790. The overshoot requirement is chosen as OS(λ, λ_(d))<20%, which lead to the parameter setting λ=3.7413 and λ_(d)=8.0097. To consider the possible change of working conditions, a set point change of 2 lbs/1000 ft² is made at t=0 s, an input disturbance of 15.5 gpm is introduced at t=1500 s, an output disturbance of 2 lbs/1000 ft² is in effect at t=3000 s, and the measurement noise is taken as zero-mean Gaussian noise with a variance of 0.414. This is shown in FIG. 12, where a graph 1200 plots the conditioned weight of a web over time (line 1202) and a graph 1204 plots the stock flow during manufacture of the web over time (line 1206). A line 1208 denotes the setpoint for the conditioned weight. From FIG. 12, it is shown that despite the large model mismatch and measurement noise, the system output robustly tracks the setpoint for all changes of working conditions, which further indicates the efficiency of the obtained parameter tuning algorithm.

Expressions of ξ, K_(x), K_(u-1), K_(xt), and K_(yt)

$\begin{matrix} {{\xi = {\left\lbrack {1,0,\ldots \mspace{14mu},0} \right\rbrack \in {\mathbb{R}}^{1 \times H_{u}}}},{K_{x} = {{- H^{- 1}}P_{cu}^{T}Q_{1}P_{cx}}},{K_{u - 1} = {{- H^{- 1}}P_{cu}^{T}Q_{3}S_{p}}},{K_{yt} = {H^{- 1}P_{cu}^{T}Q_{1}}},{K_{xt} = {H^{- 1}S_{f}^{T}Q_{3}}},{where}} & (26) \\ {{H = {{P_{cu}^{T}Q_{1}P_{cu}} + Q_{2} + {S_{f}^{T}Q_{3}S_{f}}}},} & (27) \\ {{S_{p} = {\begin{bmatrix} 1 & 1 & \ldots & 1 \end{bmatrix}^{T} \in {\mathbb{R}}^{H_{u} \times 1}}},} & (28) \\ {{S_{f} = {\begin{bmatrix} 1 & 0 & \ldots & 0 \\ 1 & 1 & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 1 & \ldots & 1 & 1 \end{bmatrix} \in {\mathbb{R}}^{H_{u} \times H_{u}}}},} & (29) \\ {{P_{cu} = \begin{bmatrix} {C_{a}B_{a}} & 0 & \ldots & 0 \\ {C_{a}A_{a}B_{a}} & {C_{a}B_{a}} & \ddots & \vdots \\ \vdots & \vdots & \ddots & \vdots \\ {C_{a}A_{a}^{H_{p} - 1}B_{a}} & {C_{a}A_{a}^{H_{p} - 2}B_{a}} & \ldots & {C_{a}A_{a}^{H_{p} - H_{u}}B_{a}} \end{bmatrix}},} & (30) \\ {P_{cx} = {\begin{bmatrix} {A_{a}^{T}C_{a}^{T}} & {\left( A_{a}^{T} \right)^{2}C_{a}^{T}} & \ldots & {\left( A_{a}^{T} \right)^{H_{p}}C_{a}^{T}} \end{bmatrix}^{T}.}} & (31) \end{matrix}$

Construction of W(S)

By using Padé approximation, Equation (18) can be changed into a real-rational form such as:

$\begin{matrix} {{\left( {P - G_{0}} \right)G_{0}^{- 1}} \approx \frac{{{g\left( {1 + {\frac{t_{d\; 0} - t_{d}}{2}s}} \right)}\left( {{T_{p\; 0}s} + 1} \right)} - {{g_{0}\left( {{T_{p}s} + 1} \right)}\left( {1 - {\frac{t_{d\; 0} - t_{d}}{2}s}} \right)}}{{g_{0}\left( {{T_{p}s} + 1} \right)}\left( {1 - {\frac{t_{d\; 0} - t_{d}}{2}s}} \right)}} & (32) \end{matrix}$

An upper bound of |(P−G₀)G₀ ⁻¹| can be found with inspection of the denominator and the numerator of Equation (32). For the numerator, the following can be calculated:

$\begin{matrix} {{{{{g\left( {1 + {\frac{t_{d\; 0} - t_{d}}{2}s}} \right)}\left( {{T_{p\; 0}s} + 1} \right)} - {{g_{0}\left( {{T_{p}s} + 1} \right)}\left( {1 - {\frac{t_{d\; 0} - t_{d}}{2}s}} \right)}}}^{2} = {{\left\lbrack {\left( {g - g_{0}} \right) + {\left( {{gT}_{p\; 0} + {g_{0}T_{p}}} \right)\frac{t_{d} - t_{d\; 0}}{2}\omega^{2}}} \right\rbrack^{2} + {\left\lbrack \left( {{gT}_{p\; 0} + {g_{0}T_{p}} + {\left( {g + g_{0}} \right)\frac{t_{d\; 0} - t_{d}}{2}}} \right) \right\rbrack^{2}\omega^{2}}} \leq {\left\{ {{\Delta \; g} + {\left\lbrack {{\left( {g_{0} + {\Delta \; g}} \right)T_{p\; 0}} + {g_{0}\left( {T_{p\; 0} + {\Delta \; T_{p}}} \right)}} \right\rbrack \omega^{2}}} \right\}^{2} + {\left\lbrack {{\Delta \; {g\left( {T_{p\; 0} + {0.5\Delta \; t_{d}}} \right)}} + {g_{0}\left( {{\Delta \; T_{p}} + {\Delta \; t_{d}}} \right)}} \right\rbrack^{2}{\omega^{2}.}}}}} & (33) \end{matrix}$

For the denominator, the following can be calculated:

$\begin{matrix} {{{{g_{0}\left( {{T_{p}s} + 1} \right)}\left( {1 - {\frac{t_{d\; 0} - t_{d}}{2}s}} \right)}}^{2} = {{{g_{0}^{2}\left\lbrack {1 + \left( {T_{p}\omega} \right)^{2}} \right\rbrack}\left\lbrack {1 + \left( {\frac{t_{d\; 0} - t_{d}}{2}\omega} \right)^{2}} \right\rbrack} \geq {g_{0}^{2}\left\lbrack {1 + {\left( {T_{p\; 0} - {\Delta \; T_{p}}} \right)^{2}\omega^{2}}} \right\rbrack}}} & (34) \\ {\mspace{79mu} {{W^{*}(s)} = \frac{\begin{matrix} {{\Delta \; g} + {\left\lbrack {{\Delta \; {g\left( {T_{p\; 0} + {0.5\Delta \; t_{d}}} \right)}} + {g_{0}\left( {{\Delta \; T_{p}} + {\Delta \; t_{d}}} \right)}} \right\rbrack s} -} \\ {\left\lbrack {{\left( {g_{0} + {\Delta \; g}} \right)T_{p\; 0}} + {g_{0}\left( {T_{p\; 0} + {\Delta \; T_{p}}} \right)}} \right\rbrack 0.5\Delta \; t_{d}s^{2}} \end{matrix}}{g_{0}\left\lbrack {{\left( {T_{p\; 0} - {\Delta \; T_{p}}} \right)s} + 1} \right\rbrack}}} & (31) \end{matrix}$

Thus, an expression of W(s) can be constructed as the formula of:

$\begin{matrix} {{W^{*}(s)} = \frac{\begin{matrix} {{\Delta \; g} + {\left\lbrack {{\Delta \; {g\left( {T_{p\; 0} + {0.5\Delta \; t_{d}}} \right)}} + {g_{0}\left( {{\Delta \; T_{p}} + {\Delta \; t_{d}}} \right)}} \right\rbrack s} -} \\ {\left\lbrack {{\left( {g_{0} + {\Delta \; g}} \right)T_{p\; 0}} + {g_{0}\left( {T_{p\; 0} + {\Delta \; T_{p}}} \right)}} \right\rbrack 0.5\Delta \; t_{d}s^{2}} \end{matrix}}{g_{0}\left\lbrack {{\left( {T_{p\; 0} - {\Delta \; T_{p}}} \right)s} + 1} \right\rbrack}} & (31) \end{matrix}$

which satisfies |W*|≧|(P−G₀)G₀ ⁻¹|. However, it can be observed that W*(s) shown above is a bit conservative since it has large values in the high frequency domain. To fix this issue, the following can be used:

$\begin{matrix} {{{W(s)} = {\frac{1}{1 - {as}}{W^{*}(s)}}},} & (32) \end{matrix}$

which results in:

$\begin{matrix} {{{W(\infty)}} = {{\frac{\left\lbrack {{\left( {g_{0} + {\Delta \; g}} \right)T_{p\; 0}} + {g_{0}\left( {T_{p\; 0} + {\Delta \; T_{p}}} \right)}} \right\rbrack 0.5\Delta \; t_{d}}{{ag}_{0}\left( {T_{p\; 0} - {\Delta \; T_{p}}} \right)}}.}} & (33) \end{matrix}$

Moreover, it can be seen that:

$\begin{matrix} {{{\left( {{P(\infty)} - {G_{0}(\infty)}} \right){G_{0}^{- 1}(\infty)}}} \approx {\frac{{gT}_{p\; 0} + {g_{0}T_{p}}}{g_{0}T_{p}}} \leq {\frac{\left\lbrack {{\left( {g_{0} + {\Delta \; g}} \right)T_{p\; 0}} + {g_{0}\left( {T_{p\; 0} + {\Delta \; T_{p}}} \right)}} \right\rbrack}{g_{0}\left( {T_{p\; 0} - {\Delta \; T_{p}}} \right)}}} & (34) \end{matrix}$

To have |W(∞)|≧|(P(∞)−G₀(∞))G₀ ⁻¹(∞)|, let a=0.5Δt_(d), and therefore:

$\begin{matrix} {{W(s)} = {\frac{1}{1 - {0.5\Delta \; t_{d}s}}{{W^{*}(s)}.}}} & (35) \end{matrix}$

CONCLUSION

In this patent document, systematic procedures have been introduced for parameter tuning of model-based controllers. Among other things, the algorithms can be efficiently employed to find a pair of λ and λ_(d) values that satisfy a user-defined specification on overshoot while guaranteeing a satisfactory settling time. The algorithms can also be extended to meet specifications on settling time. In this regard, the algorithms can ensure that a specification on settling time is achievable by a certain combination of λ and λ_(d) values, which could involve exploring the reachability set that corresponds to all stabilizing tuning parameter combinations subject to user-specified process uncertainty levels.

Note that while the techniques and algorithms described above were made with reference to tuning an MD-MPC controller, the same or similar techniques could be used to tune any other suitable model-based controller. Also, the same or similar techniques could be used to tune any suitable model-based controller in any suitable industry, not merely the paper manufacturing industry.

In some embodiments, various functions described above are implemented or supported by a computer program that is formed from computer readable program code and that is embodied in a computer readable medium. The phrase “computer readable program code” includes any type of computer code, including source code, object code, and executable code. The phrase “computer readable medium” includes any type of medium capable of being accessed by a computer, such as read only memory (ROM), random access memory (RAM), a hard disk drive, a compact disc (CD), a digital video disc (DVD), or any other type of memory. A “non-transitory” computer readable medium excludes wired, wireless, optical, or other communication links that transport transitory electrical or other signals. A non-transitory computer readable medium includes media where data can be permanently stored and media where data can be stored and later overwritten, such as a rewritable optical disc or an erasable memory device.

It may be advantageous to set forth definitions of certain words and phrases used throughout this patent document. The terms “application” and “program” refer to one or more computer programs, software components, sets of instructions, procedures, functions, objects, classes, instances, related data, or a portion thereof adapted for implementation in a suitable computer code (including source code, object code, or executable code). The terms “transmit” and “receive,” as well as derivatives thereof, encompass both direct and indirect communication. The terms “include” and “comprise,” as well as derivatives thereof, mean inclusion without limitation. The term “or” is inclusive, meaning and/or. The phrase “associated with,” as well as derivatives thereof, may mean to include, be included within, interconnect with, contain, be contained within, connect to or with, couple to or with, be communicable with, cooperate with, interleave, juxtapose, be proximate to, be bound to or with, have, have a property of, have a relationship to or with, or the like. The phrase “at least one of,” when used with a list of items, means that different combinations of one or more of the listed items may be used, and only one item in the list may be needed. For example, “at least one of: A, B, and C” includes any of the following combinations: A, B, C, A and B, A and C, B and C, and A and B and C.

While this disclosure has described certain embodiments and generally associated methods, alterations and permutations of these embodiments and methods will be apparent to those skilled in the art. Accordingly, the above description of example embodiments does not define or constrain this disclosure. Other changes, substitutions, and alterations are also possible without departing from the spirit and scope of this disclosure, as defined by the following claims. 

What is claimed is:
 1. A method comprising: obtaining information identifying uncertainties associated with multiple parameters of a model for an industrial model-based controller; obtaining information identifying multiple tuning parameters for the controller; and generating a graphical display identifying (i) one or more expected step responses of an industrial process that are based on the tuning parameters of the controller and (ii) an envelope around the one or more expected step responses that is based on the uncertainties associated with the parameters of the model.
 2. The method of claim 1, wherein the parameters comprise a process gain, a time constant, and a time delay associated with the model.
 3. The method of claim 2, wherein the uncertainties associated with the parameters of the model comprise, for each parameter of the model, an uncertainty expressed in the time domain.
 4. The method of claim 1, wherein obtaining the information identifying the tuning parameters comprises obtaining a settling time and an overshoot associated with the controller.
 5. The method of claim 4, further comprising: determining the one or more expected step responses based on a reference tracking performance ratio and a disturbance rejecting performance ratio associated with the controller.
 6. The method of claim 5, wherein determining the one or more expected step responses comprises using an iterative process to identify a value of the reference tracking performance ratio and a value of the disturbance rejecting performance ratio using a constrained optimization problem.
 7. The method of claim 6, wherein the iterative process comprises: tuning the disturbance rejecting performance ratio to obtain an optimal settling time for a fixed value of the reference tracking performance ratio; tuning the reference tracking performance ratio to activate a constraint based on the overshoot; and optimizing the settling time with respect to the reference tracking performance ratio.
 8. The method of claim 1, further comprising: identifying, for different combinations of values for the parameters of the model, multiple expected step responses of the industrial process at multiple time instances; and identifying the envelope by identifying, for each time instance, a maximum value of the step responses at that time instance and a minimum value of the step responses at that time instance; wherein the envelope defines a range of step responses given the uncertainties associated with the parameters of the model.
 9. The method of claim 1, wherein obtaining the information identifying the uncertainties associated with the parameters of the model comprises receiving the information identifying the uncertainties from a user.
 10. An apparatus comprising: at least one memory configured to store (i) information identifying uncertainties associated with multiple parameters of a model for an industrial model-based controller and (ii) information identifying multiple tuning parameters for the controller; and at least one processing device configured to generate a graphical display that identifies (i) one or more expected step responses of an industrial process that are based on the tuning parameters of the controller and (ii) an envelope around the one or more expected step responses that is based on the uncertainties associated with the parameters of the model.
 11. The apparatus of claim 10, wherein: the parameters comprise a process gain, a time constant, and a time delay associated with the model; and the uncertainties associated with the parameters of the model comprise, for each parameter of the model, an uncertainty expressed in the time domain.
 12. The apparatus of claim 10, wherein: the information identifying the tuning parameters comprises a settling time and an overshoot associated with the controller; and the at least one processing device is configured to determine the one or more expected step responses based on a reference tracking performance ratio and a disturbance rejecting performance ratio associated with the controller.
 13. The apparatus of claim 12, wherein the at least one processing device is configured to determine the one or more expected step responses using an iterative process to identify a value of the reference tracking performance ratio and a value of the disturbance rejecting performance ratio using a constrained optimization problem.
 14. The apparatus of claim 13, wherein the at least one processing device is configured during the iterative process to: tune the disturbance rejecting performance ratio to obtain an optimal settling time for a fixed value of the reference tracking performance ratio; tune the reference tracking performance ratio to activate a constraint based on the overshoot; and optimize the settling time with respect to the reference tracking performance ratio.
 15. The apparatus of claim 10, wherein the at least one processing device is configured to: identify, for different combinations of values for the parameters of the model, multiple expected step responses of the industrial process at multiple time instances; and identify the envelope by identifying, for each time instance, a maximum value of the step responses at that time instance and a minimum value of the step responses at that time instance; wherein the envelope defines a range of step responses given the uncertainties associated with the parameters of the model.
 16. A method comprising: obtaining information identifying uncertainties associated with multiple parameters of a model for an industrial model-based controller; and identifying multiple tuning parameters for the controller based on the uncertainties and one or more tuning specifications, wherein the one or more tuning specifications include constraints on a settling time and an overshoot associated with the controller.
 17. The method of claim 16, wherein identifying the multiple tuning parameters comprises identifying the multiple tuning parameters based on a reference tracking performance ratio and a disturbance rejecting performance ratio associated with the controller.
 18. The method of claim 17, wherein identifying the multiple tuning parameters further comprises using an iterative process to identify a value of the reference tracking performance ratio and a value of the disturbance rejecting performance ratio using a constrained optimization problem.
 19. The method of claim 18, wherein the iterative process comprises: tuning the disturbance rejecting performance ratio to obtain an optimal settling time for a fixed value of the reference tracking performance ratio; tuning the reference tracking performance ratio to activate the constraint on the overshoot; and optimizing the settling time with respect to the reference tracking performance ratio.
 20. The method of claim 16, wherein: the parameters comprise a process gain, a time constant, and a time delay associated with the model; and the uncertainties associated with the parameters of the model comprise, for each parameter of the model, an uncertainty expressed in the time domain. 